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ABSTRACT 

We examine the effect of different radiative transfer schemes on the properties of 3D simulations of 
near-surface stellar convection in the superadiabatic layer, where energy transport transitions from 
fully convective to fully radiative. We employ two radiative transfer schemes that fundamentally differ 
in the way they cover the 3D domain. The first solver approximates domain coverage with moments, 
while the second solver samples the 3D domain with ray integrations. By comparing simulations that 
differ only in their respective radiative transfer methods, we are able to isolate the effect that radiative 
efficiency has on the structure of the superadiabatic layer. We find the simulations to be in good 
general agreement, but they show distinct differences in the thermal structure in the superadiabatic 
layer and atmosphere. 
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1. INTRODUCTION 

Techniques for constructing ID stellar models continu- 
ally evolve in order to encompass more realistic physics. 
A stellar structure model is computed with a given set 
of physical assumptions and boundary conditions. Stel- 
lar evolution models with convective envelopes require 
an approximation for convection to define the structure 
within the convection zone. Current stellar models with 
convective envelopes are not accurate near the surface. 
This is largely due to the lack of a proper description for 
inefficient convection. 

Th e Mixing Length Theory (MLT; iBohm-Vitensel 
1958) is the most commonly used approximation for con- 
vection in ID stellar models. MLT represents the range 
of turbulent eddy-sizes with a single scale length that is 
proportional to the pressure scale height. The so-called 
"mixing length" is I = aHp, where a is a free parameter. 
This prescription is able to accurately model the stellar 
structure where convection is efficient and the tempera- 
ture gradient is nearly adiabatic. In the superadiabatic 
layer (SAL) near the surface, however, the energy trans- 
port transitions from convective to radiative and MLT 
breaks down . Note that some MLT treatments (see e.g., 
lArnett et al.ll2010h include an additional free parameter 
that sets the ratio between the eddy size and the dis- 
tance it travels. Other appr oximations for convec t ion ex - 
ist, such as the treatment of lCanuto fc Mazzitellil (|1991l ). 
but it, too, contains a length scale which has to be mod- 
eled. In this region of inefficient convection and convec- 
tive overshoot we are unable to accurately model stellar 
structure. 

One of the major limitations of using MLT or similar 
treatments of convection in stellar models is that it effec- 
tively renders the stellar radius to be a free parameter. In 
ID models with a given set of physics, the mixing length 
parameter a determines the specific entropy of the con- 
vection zone, which in turn sets the radius of the model. 
In addition to the approximation for convection, stel- 
lar models require an atmospheric boundary condition. 
Many stellar evolution codes force the atmospheric struc- 
ture to be the Eddington T-t relation, although other 
model atmospheres can be used as well. In the case of 
the Eddington T-t relation, the atmosphere is purely 



radiative, and does not include the effect of convective 
overshoot. This, too, contributes to the inaccuracy of 
near-surface layers in stellar models. 

A promising strategy for improving our understand- 
ing of inefficient convection is to use 3D radiation hy- 
drodynamic (RHD) large eddy simulations (LES) to 
calculate a more realist ic stratification. Pioneered by 
iNordlundl (|1982t Il985b| ). such simulations account for 
both the realistic transfer of energy through convec- 
tive gas dynamics and radiative transfer in the SAL 
where energy transport changes from convection to ra- 
diation. Simulations provide a realistic stratification 
that sclf-consistently couples the convective envelope 
to the radiative atmosphere. Simulating convection 
has been fe asible for several decades, and a number of 
groups le.g iChan fe Sofialll989HStein fc Nordlvmdlll989L 
19981 120001: iCattaneo et all 119911: iPorter fc Woodward! 
1994 iKim et alJll995t iRobinson et al II2003L 12004 120051: 
Jung et al.l 120071 iFrevtag et al.1 l20TTl etc.) have per- 
formed 3D simulations of convection. Simulations 
have been used to study 3D penetrative convec- 
tion both above and below convectiv e envelopes (e.g. 



Singh. Roxburgh fc Chanl 11994 119951: iMuthsam etafl 
1993) and have also been compared with ID stellar 
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and atmos phere models (e.g. INordlund fc Dravins 
Kim et al.lll996t IFrevtag et al.lll996t lAbbett et aL._ 
Nordlund fc Steinl 120001: IPiau fc Steinl 120081: iPiau et al.l 
2011aHbl ld)~ 

Several codes are currently being used by various 
groups to simulate near surf ace convection. These 
codes include STAGGER (INordlund fc Galsgaard" 
[1995[ IGalsgaard fc Nordlundl 119961: IStein fc Nordlunc 
1 1998 1 : lHavek et al l I2011D, CQ5BO LD (IFrevtag et ~ 
120021 l20ll . MURaM lVogler et al J 120051). ANTARES 
(IMuthsam et al.l 120071 120101), and our code, which is 
based on the code o f lKimfc Char](fT998h . Each of these 
codes has been developed independently, and can be 
quite different in their numerical treatment of radiative 
transfer, subgrid scale models, advection schemes, as 
well as input physics such as equations of state and 
opacity. 

Simulations have been successfully applied to the outer 
convection zone in the Sun and in a select group of in- 
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dividual stars. Such simulations have led to improved 
p-mode oscillation frequencies in the outer layers of the 
models where errors in the sound speed contribute to dis- 
crepancies between the observed and model frequencies. 
Including some turbulent effects has l ead to improved 
frequencies in stan dard solar models ([Rosenthal et al.1 
1999 : ILi et al ll2002l) as w ell as rj Bootis and a Centauri 
( Straka et al.ll2006L [2007^ . 

Simulating individual stars has been useful for under- 
standing how convection affects stellar structure in par- 
ticular cases, and a natural extension of this strategy is 
to simulate surface convection across the HR diagram. 
Some progress has been made in parameterizing the tur- 
bule nt dynamics i n the SAL in thi s wa y. For exam- 
ple. iLudwig et al.l (|1995l 119981 I1999D and iFrevtag et all 
(1999) have undertaken extensive efforts to map the mix- 
ing efficiency (or effective a for MLT) using 2D RHD 
simulations, and constructing ID envelope models by 
matching the specifi c entropy of the mod e ls an d sim- 
ulations. Recently, iTrampedach fc Stein! (|2011[ ) used 
3D simulations to map convective efficiency in a range 
of the HR di agram that include s main -sequence stars 
and giants. ITrampedach et all (| 19991 ) have also at- 
tempted to extract a mixing length parameter from sim- 
ulations by matching averaged 3D simulation stratifica- 
tions to ID envelope models and modifying the turbu- 
lent pressure. Simulations have also been applied to 
the formati on of spectral l ines and int erferometric ob- 
servables. lAsplund et al.1 (2000, 2009) famously used 
3D spectral synthesis to study solar line profiles for 
abundance analysis, but similar techniques have since 
been used on the sun (e.g. ICaffau et al.ll20Tol l20ll and 



other stars (e.g. iRaimrez et al.l 12009 ; Collet et aTT 



2009: iChiavassa et al.l I2009L 120101 : IWende et all 
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2009; 



Behara et al.ll2010h . 



While these applications have been largely successful, 
the accuracy of simulations is limited by what physi- 
cal processes are included, and in what manner they 
are represented numerically. One way to get a handle 
on the uncertainty is through comparing different ap- 
proximations for the same physical process. Encourag- 
ingly, comparisons by iKupkal <|2009f ) , and more recently 
iBeeck et al.l (|2012f) . show that existing codes are in good 
agreement for simulations of the sun. There are, how- 
ever, apparent differences in the thermal structure. Be- 
cause there are numerous differences between each of the 
codes, the source of the differences cann ot be tr a ced to 
a p articular method or ap proximation. iVoglerl (|2004l ) 
and IChiavassa et al.l ()201lD have looked at the effect of 
gray and non-gray opacities in 3D simulations, but their 
comparisons necessarily used the same radiative transfer 
method. 

The goal of this work is to examine the extent to which 
the radiative transfer method affects the resulting ther- 
mal structure and gas dynamics of the simulation, while 
keeping the rest of the physics the same. In Sections [5] 
and [3] we describe our code and the two radiative trans- 
fer schemes that we have implemented. Section U com- 
pares results from simulations computed with each of the 
solvers. 

2. THE RADIATION HYDRODYNAMICS CODE 

Stellar convection zones can span many scale heights 
(the solar convection zone exceeds 20 pressure scale 



heights) and evolve over a range in timescales. The ther- 
mal relaxation time of an entire convection zone is too 
large to be simulated entirely with local RHD LES, so 
the simulation domain is limited to several scale heights 
near the surface. Energy transport transitions from fully 
convective at the bottom of the computational domain, 
to fully radiative at the top. 

The radiation h ydrod ynamics code is based on the code 
of iKim fc Chanl (|19987). It solves the equations o f hy- 
drodynamics as described bv lRobinson et all (|2003l ). but 
the radiative term Q ra d is computed using either the 3D 
Eddington approximation or long-characteristic ray inte- 
gration. Radiative transfer is treated with the diffusion 
approximation in the lower part of the domain (r > 10 4 ) 
where it is valid. 

Two popular simplifications to the governing equations 
are the Boussinesq and anelastic approximations. In 
both cases, sound waves are filtered out which relaxes 
the time step restriction imposed by the CFL condition. 
Neither of these two approximations are suitable for sim- 
ulating realistic stellar surface convection. The Boussi- 
nesq approximation is not valid because it does not per- 
mit strong stratification and the Mach number must be 
small. The anelastic approximation, which can be used 
in a heavily stratified flows, is valid only in the deep con- 
vection zone but is unsuitable for convection near the 
SAL where the Mach number is high. Instead, we solve 
the fully compressible Navier-Stokes equations: 
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where E — e + pv 2 /2 is the total energy density and p, 
v, P, e, and g are the density, velocity, pressure, specific 
internal energy, and gravitational acceleration, respec- 
tively. The viscous stress tensor for a Newtonian fluid is 
Sij = p{dvi/dxj + dvj/dxi) — 2/3/i(V • v)Sij. The radia- 
tive term Q ra d is computed using one of the two radiative 
transfer schemes described in Section [3] Explicitly cou- 
pling the radiative term, Q ra d, as a source term in the 
hydrodynamics equations imposes a restriction on the 
time step. Switching radiative transfer solvers changes 
the method for computing Q ra d, but does not ease the 
restriction. 

The nature of stellar convection makes it impossi- 
ble to directly simulate the range in dynamical scales. 
LES resolve the scales where most of the turbulent en- 
ergy transport takes place, however, it is necessary to 
account for the energy transport on unresolved scales. 
We accomplish this by adjusting the dynamic viscos- 
ity to account for sub- grid scale stresses according to 
ISmagorinskvl (|1963D . 

Where possible, we keep the physics of the simula- 
tion con sistent with that of the ID stellar evolution code 
YREC (jDemarque et al.l [2008h . T he code uses OPAL 
equat ion of sta te and opacity tables (iRogers fc Navfonovl 
I2002T) . and the lFerguson et al.l 1)20051 ) tables for low tem- 
perature opacities. The two simulations presented in 
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this work use Rosseland mean opacity computed from 
the tables. The opacit y and equation of state a re both 
for the GS98 mixture (|Grevesse fc Sauvall [l998h . Fur- 
ther, we treat radiative transfer in local thermodynamic 
equilibrium (LTE) with gray opacities, but these are not 
limitations of the methods; both can be extended to ac- 
count for non-g r ay effects using; met hods described by 
INordlundl (|1982f ). lLudwig et all <[199l . and lVogler et al.l 
(120041) . 

3. RADIATIVE TRANSFER IN 3D SIMULATIONS 

There are numerous methods for t he treatment of mul- 
tidimensional radiative transfer (see lCarls son (2008]) and 
references therein for a review of various met hods). Ra- 
diative transfer in 3D was first presented by INordlundl 
(|1985aD . Today, the most com monly used techniques for 
3D radiative transfer arc long- (ICa^oiJl^lAuerlT96l 
and short- (jKunasz fc Auerl Il988f ) characteristic meth- 
ods. Both of these methods share the same strategy for 
3D spatial coverage, wherein they combine ID solutions 
along a few selected rays to sample all An steradians. In 
the following sections we describe the methods of long- 
characteristic ray integration and the 3D Eddington ap- 
proximation, which represent two fundamentally differ- 
ent approaches to achieve spatial coverage. Domain cov- 
erage is important because the transfer of radiation is 
not a local problem in the SAL. 

Perhaps the most significant challenge for any 3D ra- 
diative transfer solver is to adequately account for the lo- 
cal angular variation of Qrad (see Section 13. 2. 2 [) . The two 
methods that we compare approximate 3D domain cov- 
erage in very different ways. Ray integration techniques 
estimate <5 ra d by solving the transfer equation along ID 
characteristics, but they have potentially poor spatial 
coverage that is achieved through combining a few ID 
solutions. Achieving complete domain coverage in this 
manner is not computationally feasible, so characteris- 
tic methods are necessarily undersampled. The 3D Ed- 
dington approximation approaches the problem entirely 
differently. In it, radiative transfer and domain coverage 
are approximated at the outset through moment equa- 
tions that are closed at the first order. In the follow- 
ing sections we investigate the extent to which these two 
approximations affect the structure of the SAL. The ac- 
curacy of the spatial coverage in moment methods such 
as the 3D Eddington approximation can be enhanced by 
taking higher order moments, which is analogous to im- 
proving the spatial coverage of ray integration methods 
with additional characteristics. 

3.1. 3D Eddington Approximation 

As an alternative to characteristic based ray integra- 
tion methods, the radiative transfer equation can be ap- 
proximated through m oment equations as described by 
lUnno fc Spiegel! (fl966h : 

V- | — VJ) -pnJ + P KB = 0. (4) 
\3pn J 

and the net radiative heating or cooling is: 

Qrad = 4 P K(J - B) (5) 

The 3D Eddington approximation has been us ed in 
codes applied to other astrophysical regimes (e.g. iBossI 



H980l Ii98l . although to our knowledge our code is 
the only application of it to simulations of near-surface 
convection. The Eddington approximation describes 
isotropic radiation exactly in both optically thin and 
thick regions, and is not equivalent to the diffusio n 
approximation in optically thin regions (|Ruttenll200l . 
Further, although we use the Planck function as the 
source function, the 3D Eddington approximation does 
not strictly require LTE. 

The 3D Eddington approximation has advantages over 
other methods, afforded by its simplicity and speed. One 
of the advantages of the 3D Eddington approximation 
is that the nature of its formulation alleviates concern 
that a nearby source or sink of energy will be missed 
as a result of limited domain coverage. This is not to 
say that the spatial domain is fully sampled (it is still 
defined by first order moments), but it is complete in 
that it includes information from all other points and 
not just those that would be sampled by characteristics. 
The angular resolution is limited by the first order mo- 
ments, but it is sensitive to the complete domain be- 
cause of the iterative nature of the numerical solvers. 
At a particular location at an instant in time, the lo- 
cal (5 ra d computed from ray integration only includes 
information from along the characteristics, which could 
potentially miss an important source or sink. Two grid 
points that are not connected by a ray can only influ- 
ence each other indirectly. In simulations of convection 
in the SAL, this risk is largely mitigated by the time evo- 
lution, or by rotatin g the undersampled bund le of rays in 
a manner similar to lStein fc Nordlundl (|2003[) . Complete 
domain sensitivity can be ensured by iterating on the ray 
solution, but this adds to the already high computational 
cost. The 3D Eddington approximation achieves domain 
coverage without requiring such methods. 

The second advantage is that the 3D Eddington ap- 
proximation can be solved with a simple and fast numer- 
ical method. This is in a large part because the solution is 
computed directly on the hydrodynamic grid, and there 
is no need to define additional grids for radiative transfer 
along each characteristic. Aside from the computational 
convenience of having the radiation and hydrodynamic 
solvers share the same grid, errors introduced from inter- 
polation are also avoided. Interpolation is likely not the 
dominant source of error, and is almost certainly smaller 
than the error introduced by the 3D Eddington approxi- 
mation for spatial coverage or the undersampled numer- 
ical quadrature, so avoiding interpolation is primarily a 
performance advantage. Interpolations in ray integra- 
tions can be avoided altogether by casting rays so that 
they always intersect the hydrodynamic grid, but flexibil- 
ity in ray placement would be lost, resulting in reduced 
domain sensitivity. 

3.1.1. Iterative Solvers 

Solving large sets of partial differential equations has 
been dealt with in great detail, and there exists a rich 
body of work from which we can develop solvers for 
the anisotropi c inhom ogeneous Helmholtz equation (see 
iDouglas et al.l (|2003[ ) and references therein). The 3D 
Eddington approximation, in conjunction with assuming 
LTE to set the source function, results in a linear rela- 
tionship for the mean intensity, J. In principle, Equation 
(j4]) can be written as a large linear system and solved di- 
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rectly, but this is not done in practice because of the 
enormous computational cost. Much faster iterative (re- 
laxation) methods are used instead. 

While classical relaxation techniques are sufficient to 
solve Equation (|4j, methods with better convergence 
properties, such as alternating-direction-implicit (ADI) 
and multigrid methods, reduce the computational bur- 
den significantly. Multigrid methods are particularly 
powerful, as they can render the comp utational cost pro - 
portional to the number of grid points ([Press et al.f l992). 
In the following sections we describe two solvers that we 
use in the code. Note that the two methods are exactly 
equivalent in that they produce solutions to Equation 
j4} that agree to within the specified error threshold. 
Because the solvers result in the same solution, select- 
ing a solver is based on numerical stability and parallel 
performance considera tions. 

First introduced bv IDouglas k, Peacemanl (|1955l ) and 
iPeaceman fc Rachfordl ( ~955l ). ADI methods are use- 
ful on regularly spac ed multidimensional grids. Both 
iPress et all (|1992t ) and IDouglas et al.l (|2003t ) describe the 
method in detail. Each ADI iteration is comprised of 
three one-dimensional implicit problems, obtained by 
splitting the terms in equation |4] according to spatial di- 
mension. For each spatial dimension, the derivatives are 
treated implicitly only for that dimension. The terms 
corresponding to each dimension represent a series of 
tridiagonal l inear systems, for whi ch there are many fast 
solvers (e.g. lAnderson et aLlll999| ). 

The systems are solved sequentially in 'sweeps' along 
each dimension. Completing the three directional sweeps 
provides an improved estimate of the solution. The final 
solution is achieved by iterating until the desired error 
threshold is reached. The convergence properties of the 
ADI technique are sensitive to the initial guess, asym- 
metry of the flow, and the problem size. Our tests show 
that convergence is not always linear with problem size, 
which can become cost-prohibitive for very large grids. 

Unlike the ray integration scheme, the ADI solver 
does not present a natural parallelism that is as easy 
to exploit. Although each tridiagonal system within a 
directional sweep is independent, the successive itera- 
tions can result in very high communication costs. To 
achieve scalability on distributed memory machines, we 
use a process-scheduling algorithm with asynchronous 
communication to overlap computation and communi- 
cation. While it is not our intent to discuss p a rallel algo- 
rithms in deta i l, we point out that iPovitskvl (|2002fl and 
IDouglas et al.l (|2003f ) discuss strategies for parallel ADI 
algorithms. 

We have also implemented a so-called optimal algo- 
rithm using a geometric multigrid method. Such meth- 
ods solve partial differential equations using a recursive 
procedure applied to a hierarchy of grids. Multigrid 
solvers achieve better convergence rates by reducing low- 
frequency errors on c oarser grids and int erpolating the 
solution to finer grids ()Douglas et al.ll2003l ). The solution 
is mapped between grids using smoothing, restriction, 
and interpolation operators. The two major advantages 
of this method are stability and scalability. 

The multigrid method is much more stable than ADI 
in regions of anisotropic flow. The ADI method requires 
a good initial guess for the solution. Using the solution 
at the previous time step for the initial guess is usually 



is sufficient, but if the flow is highly anisotropic, the ADI 
method can fail to converge. Unlike the ADI solver, the 
multigrid solver does not require a good initial guess for 
the solution. This makes the solver very stable, even for 
highly turbulent flows, such as granulation in giants and 
early type stars. 

The scalability of the solver is optimal in the sense that 
in scales linearly with problem size. That is, is solves 
the equation with operational cost proportional to the 
number of unknowns. Our tests confirm that the solution 
is achieved with O(N) operations without a good initial 
guess. Although not strictly required, a good initial guess 
improves the speed of the solver. 

Our multigrid m ethod is bas e d on the 'full multigrid' 
solver described bv IPress et al.l (|1992l ). modified to solve 
the anisotropic inhomogeneous Helmholtz equation in 
3D, with periodic boundary conditions on the horizon- 
tal walls and Neumann boundary conditions at the top. 
It employs W-cycle structure for multigrid cycles, and 
Gauss-Seidel smoothing with red-black ordering. As yet, 
our multigrid solver is not parallelized for distributed 
memory, which limits its use to shared memory machines, 
however, a full parallel version of this solve r can be imple- 
mente d following the strategy outlined bv IDouglas et al.l 
(|200l . 

3.2. Long Characteristic Ray Integration 

Ray integration methods are attractive because they 
solve the integro-diffcrcntial radiative transfer equation 
along rays (or characteristics), which provides a robust 
solution along a ID path. Extending the solution to 
3D, however, requires solving the transfer equation along 
many ray paths that cover every direction. Good domain 
coverage can be achieved when many rays are used, but 
the computational cost of this is typically too high. To 
alleviate the computational burden, only a few rays are 
used. 

In contrast to the simplicity of the 3D Eddington ap- 
proximation, long-characteristic ray integration is com- 
prised of a ray integration step and a spatial integration 
step that must be carried out at each location in the 
computational domain. The first part of the solution is 
to perform the ID radiative transfer calculations at var- 
ious angles through the domain. The second part of the 
solution is to integrate the result from the ID rays at each 
grid point over the unit sphere with numerical quadra- 
ture. The details of how these two steps are implemented 
can have an effect on the final solution. We discuss our 
implementation in Sections 13.2.11 and 13.2.21 below. 

In terms of computational performance, our 3D Ed- 
dington solvers outperform the ray integrators (see Sec- 
tion 13.41 for a detailed comparison), but ray integra- 
tion presents a natural parallelism that can be ex- 
ploited. In modern computing, scalability is often 
more important than th e inh erent cost of the algori thm. 
iHeinemann et al.l (|2006D and iMuthsam etaLl (|2010l) de- 
scribe parallel algorithms for ray integration that can be 
applied to distributed memory machines using domain 
decomposition. The good performance over a large num- 
ber of distributed memory compute nodes makes ray in- 
tegration an attractive option for modern codes. 

3.2.1. Solving the Transfer Equation 
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Following the description o f iStein fc Nordlundl (|2003l) . 
we adopt the lFeautrierl (| 19641 ) representation of the trans- 
fer equation for the numerical solution: 



where S is the source function (taken to be the Planck 
function), P = 1/2 [7+ (r) - J_(r)] and J+(r) and 7_(r) 
are the outgoing and incoming specific intensities at opti- 
cal depth r. The local heating/cooling Q ra d is computed 
by integrating in each direction: 

Q + = 1+ - S, q- = T - S, (7) 

and the contribution to the net radiative heating or cool- 
ing for the ray is 

Qrad = Q + +Q~ (8) 

We use the Hermitian differencing method of lAuerl 
(1976), which achieves 4th order accuracy while main- 
taining essentially the same cost as a 2nd order scheme. 
We found similar results for the final Q ra d when using 
the 2nd and 4th order schemes, but the sensitivity on 
the result may depend on the spatial resolution of the 
simulation, as well as the method for interpolating be- 
tween the hydrodynamic and radiation grids (discussed 
below). 

We define all ray characteristics to begin at the top of 
the box and integrate into higher optical depth. Rays are 
cast in directions defined by two angles (9, (j>). The 9 an- 
gle is measured with respect to the vertical direction, and 
<p is the azimuth angle. Each ray direction corresponds to 
a bundle of parallel rays that traverse the computational 
domain. In order to compute the formal solution along 
a given characteristic, the source function and optical 
depth must be known at each point along the ray. From 
this, the radiative heating or cooling Qrad is computed at 
each of the points where the source function and optical 
depth are known. Because the ray integration is per- 
formed along lines that are not usually coincident with 
the hydrodynamic grid (the vertical ray is an exception), 
the required thermodynamic variables must be mapped 
from the hydrodynamic grid to the radiation grid. 

The hydrodynamic grid is defined by the regular 3D 
Cartesian mesh that specifies the locations of zone- 
centered quantities. The temperature, density, and opac- 
ity are tracked on this grid. The radiation grid is defined 
by the intersection points of a parallel ray bundle with 
the hydrodynamic grid. Intersection points occur in 2D 
planes between the zone centers of the hydrodynamic 
grid, so a 2D interpolation is required at each intersec- 
tion to determine the required thermodynamic quantities 
for ray integration. Figure [1] shows a schematic with one 
ray intersecting the hydrodynamic grid at a steep and 
shallow angle. 

A ray always has a vertical component (we do not use 
periodic rays), so the final solution only requires points 
in the horizontal planes of the hydrodynamic grid. Very 
shallow rays can traverse long distances between succes- 
sive intersections of the horizontal planes of the hydro- 
dynamic grid, so we include intersections of the vertical 
planes as well to ensure that each ray is adequately sam- 
pled before performing the Feautrier integration. The 
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Fig. 1. — A 2D schematic illustration of a steep (top) and shallow 
(bottom) long ray cast through a periodic domain. Mapping be- 
tween the ray and hydrodynamic grids requires interpolating data 
from the hydrodynamic grid (blue points) to the red points. The 
solution along a characteristic (red points) is only required at inter- 
sections of horizontal planes (solid red points) but all intersection 
points are used to avoid very large steps in r for shallow rays. 

intersections with vertical planes are only included if the 
distance to the next horizontal plane is large. This re- 
striction avoids having points that are spatially very close 
together. The result is that the forward mapping (from 
the hydrodynamic to radiation grid) includes most inter- 
sections of the ray with the hydrodynamic grid, but the 
reverse mapping (from the radiation to hydrodynamic 
grid) only includes intersections of the ray with horizon- 
tal planes of the hydrodynamic grid. 

The mapping between the two grids introduces a 
smoothing to the solution twice. First, the hydrody- 
namic quantities (in this case pn and T) must be in- 
terpolated to points along the ray bundles. The sec- 
ond smoothing occurs with the reverse transformation of 
interpolating the computed Q r ad back to the hydrody- 
namic grid. The smoothing effect can be somewhat mit- 
igated by employing higher order interpolation schemes, 
such as bicubic or distance weighted parabolas, but this 
will add to the computational expense. 

3.2.2. Spatial Coverage 

After the ray integration step, the contribution to the 
radiative heating or cooling is known for each ray direc- 
tion, and the result of all rays at each grid point must 
be integrated over the unit sphere. To completely cover 
the 3D domain, rays need to be cast along every line of 
sight over Air steradians, at every point in the grid. This 
is far too computationally demanding, so we use only a 
few rays which arc integrated over An steradians using a 
quadrature rule. There a re numerous quadrature rules 
av ailable ( for ex ample, see lAbramowitz fc Stegunl ()1972f ) 
or lAntial (|1991| )) and the choice of rule may affect the 
total radiative flux in the simulation. 

For local RHD simulations we can take advantage of 
azimuthal symmetry by spacing rays at fixed intervals 
of 4>, and the integration reduces to a ID quadrature 
along the azimuth, and a second ID quadrature along 
the cosine-angle. If the rays equally spaced in azimuth, 
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the weights of rays with the same 9 are equal, and the 
unit sphere integration is reduced to a one dimensional 
calculation. Using a quadrature rule, the integration is a 
weighted sum of the radiative heating at each grid point. 
Adding to the number of rays improves the accuracy of 
the unit sphere integration, and achieving good spatial 
coverage is important for accurately calculating the local 
and total energy flux. 

Typically only a few rays are used to keep the compu- 
tational expense manageable. For example, the STAG- 
GER code uses 5 rays which are rotated by 15 degrees 
every time st ep dStein fc Nordlundl I2003T ) . or 9 rays in 
recent work (|Beeck et al.l l2012f T The solution is quite 
sensitive to integration over 9 because the vertical en- 
ergy flux is the largest contributor to the total flux. If 
the ray directions are not constrained by the hydrody- 
namic grid, the theta angle of the rays can be chosen 
to correspond to a desired quadr ature rule. In this work 
we use Gauss-Radau quadrature ( Hildcbrand 1956j) with 
one vertical ray and four rays inclined at cos6> = 1/3, 
whi ch is the same as t hat u sed in the STAGGER code 
(in IStein fe Nordlundl (|2003Q ). The vertical ray has a 
weight of 1 /4 and the four inclined rays have a combined 
weight of 3/4. 

3.3. Variation of Local Radiative Heating and Cooling 

Spatial coverage is critical for radiative transfer in sim- 
ulations of near-surface convection. A radiative transfer 
solver must have good spatial coverage because of the 
large asymmetries that convection introduces in the flow, 
and because the SAL spans the transition from optically 
thick to thin where radiative transfer becomes non-local. 
The cost of good spatial coverage can be problematic 
for 3D RHD simulations, so approximations are made to 
the angular resolution to keep the computational costs 
affordable. The 3D Eddington approximation simplifies 
spatial coverage through moments, while they ray inte- 
gration relies on numerical quadrature. In both cases the 
result can be inaccurate if the function is highly variable 
or asymmetric. Both of these properties exist for the ra- 
diative heating in the SAL and pose a challenge for any 
affordable 3D radiative transfer method. 

The stellar flux is predominantly vertical (introducing 
large asymmetry in the function) and large horizontal 
variation can occur from the high contrast of hot up- 
flowing granules and cool inter-granular lanes. The way 
that radiative transfer schemes deal with this is a large 
source of uncertainty and limits the absolute accuracy 
of the methods. For example, spatial quadrature rules 
based on different sets of polynomials will assign differ- 
ent weights to the vertical and horizontal rays, which can 
result in different estimates for the total energy fluxes 
from the same opacity and source function. This is a 
significant source of uncertainty, and one that can only 
be eliminated by using a very large number of rays. 

To show how the ray integration method is sensitive 
to spatial coverage, we compute Qrad using the long- 
characteristic ray solver with 65 rays, on a snapshot from 
the thermally relaxed 3D Eddington simulation. Note 
that the number of rays is far in excess of what is typical 
or affordable for a simulation. In the simulations, the 
result of the ray calculation is integrated over the unit 
sphere to achieve the net Q ra d, which is then applied 
to the Navier Stokes equations. In this case, however, 



we average over azimuth and present the local Qrad as 
a function of cosine-angle. Figure [2] presents several ex- 
amples of the local Q ra d corresponding to three different 
regions of surface granulation at the r = 1 surface. We 
present the local <3 ra d in dimensionless units, normalized 
to the total stellar flux. Even with 65 rays, the variation 
in local Q r ad is not always well resolved. 

The curve in the center panel of Figurc[2]shows Q ra d for 
the center of a granule, where the local flux is dominated 
by the vertical ray. This is expected, as it is the hot 
upflows that contribute the majority of the stellar flux. 
Because the center of a granule is hot and surrounded by 
hot gas, the heating from horizontal directions is quite 
small relative to the vertical flux. In upflowing regions 
such as this, the large variation of the local Q r ad function 
can make it difficult to determine the local energy flux if 
only a few rays are used. 

In the narrow intergranular lanes, the gas is cool and 
surrounded by much hotter upflows. Here, the local en- 
ergy flux is still dominated by the vertical direction, but 
the magnitude of the flux is much smaller than in the 
hot upflow. Although there is still a net heat flux verti- 
cally, nearby upflows are heating the gas from shallower 
(horizontal) angles. The intergranular lanes show a large 
asymmetry of the local Q ra d function, and the local flux 
can change sign at small angles. Both of these properties 
can be troublesome for a sparsely sampled quadrature 
rule. 

Regions of mixed flow, such as a small upflow that 
is near to an intergranular lane, show characteristics of 
both the hot and cold extremes described previously. 
The mixed flow example presented in Figure [5] shows a 
similar vertical flux as in the downflow region, but has 
significant heating from angles far from vertical. Fur- 
thermore, the variation in heating as a function of angle 
varies considerably over small angular increments. The 
local flux is much smaller than that from the center of 
a large granule, but the local Q rac i function is still very 
asymmetric and dominated by the vertical direction. 

The three examples in Figure [2] show the angular vari- 
ation in Q ra d at an instant in time, and do not represent 
the bulk flow. Because of the highly variable nature of 
the local Q ra d profile, placing rays to adequately sample 
the vertical flux may be important for accurately deter- 
mining the radiative energy flux. 

3.4. Performance Considerations 

Our comparison in Section [4] of the moment-based 3D 
Eddington approximation with the under-sampled ray 
integration method shows that they result in largely 
equivalent stratifications through the stellar SAL. Cur- 
rent 3D RHD simulations are limited by computational 
power, so it is prudent to compare the costs of these 
two methods that produce equivalent results. Numeri- 
cal solvers for moment based methods are fast, and offer 
significant computational savings. 

In our experience, the 3D Eddington solvers typically 
out-perform the ray integrators by a significant margin. 
Figure [3] shows the time-to-solution for simulation do- 
mains of varying size. The radiation solvers are applied 
only to the top of the simulation domain where the tem- 
perature gradient is significantly superadiabatic. The 
cases presented here correspond to a simulations with 
a superadiabatic domains that are 65 zones deep and 
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Fig. 2. — The azimuth-averaged angular variation in the local radiative heating and cooling Q ra d (in dimensionless units normalized to 
the stellar flux) for three locations (marked in the left panel) near the optical surface. The area enclosed by each Q ra( j curve is the local 
net heating or cooling. The left panel shows the locations (large upflow in red, small upflow in green, and downflow in blue) for the local 
Qrad profiles in the other two panels. The Q ra( j profiles can exhibit large asymmetry and variation over small changes in angle. Using 
quadrature rules with two or three angles can systematically under- or over-estimate the stellar flux. 

Similarly, the cost of both the ADI and multigrid 
solvers for the 3D Eddington equation show a roughly 
linear proportionality with grid size. A good initial guess 
is necessary for the ADI method to achieve linear scaling, 
but no such guess is necessary for the multigrid method. 
The iterative 3D Eddington solvers have a tolerance pa- 
rameter that defines how close the numerical solution 
matches equation ((4]) We see a performance benefit by 
adjusting the tolerance parameter from 10~ 8 to 1CP 6 . 

The speed difference between the moment and charac- 
teristics methods can make a tangible difference. Given 
that the computational cost of a 3D radiation hydrody- 
namics simulation is typically dominated by the radiation 
solver, the transport of radiation is simplified, for exam- 
ple, by limiting the domain coverage, or by limiting the 
frequency coverage to a few opacity bins. The computa- 
tional cost savings afforded by moment based methods, 
which appear to offer similar domain coverage as under- 
sampled ray integration, could be used to better cover 
the frequency domain using opacity binning or sampling 
techniques. 
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Fig. 3. — Performance comparison of the ray integrator with the 
ADI and multigrid solvers for the 3D Eddington equation. All tests 
were performed on the same hardware using a single processor. The 
3D Eddington solvers have a clear performance advantage over ray 
integration. 



with varying horizontal sizes. All simulation runs for the 
comparison were performed on the same hardware using 
a single processor. 

The ray integration solver presented in figure [3] is our 
most economical one, which uses 5 rays with bilinear 
interpolation. Note that the cost of additional rays is 
not simply proportional to the number of rays because 
shallow rays traverse a longer path through the simula- 
tion domain. A significant improvement in performance 
could be reali z ed by adopting the scheme described by 
iFrevtag et al.1 (j2012j ). where only ray intersections with 
horizontal planes are included, and by aligning the az- 
imuthal angle of the rays with the hydrodynamic grid, 
thereby reducing interpolations to one dimension. The 
smaller computational cost of the one dimensional inter- 
polations comes with the consequence of reduced domain 
sensitivity because the azimuth angle of the rays could 
no longer be rotated each time step. Because the ray in- 
tegrator is not iterative, for a given ray configuration the 
amount of computational work is linearly proportional 
to the size of the grid. Increasing the horizontal extent 
of the simulation domain results in the expected linear 
increase in computational cost. 



4. COMPARING RADIATIVE TRANSFER SCHEMES IN 3D 
SIMULATIONS 

We have implemented the two methods described in 
Section [3] for computing radiative heating and cooling in 
a 3D LES. While differences in the methods are clear, the 
only way to determine their effect on the thermal struc- 
ture and gas dynamics from solving the Navier-Stokcs 
equations is to use them in simulations. Because Q ra d is 
essential in determining the thermal structure, we can- 
not compare methods on the same simulation snapshot 
at an instant in time. Instead, we compare statistics 
gathered from long runs of simulations computed with 
each radiative transfer solver. 

In this section we compare results of two simulations: 
one computed with the 3D Eddington approximation 
using the ADI solver (described in Section I3.1[) . and 
the other with long-characteristic ray integration (de- 
scribed in Section I3.2[) . The two simulations have the 
same composition (Z = 0.0169, X = 0.7366), energy 
flux (T cff « 5725K), and surface gravity (log(g) = 4.438). 
Both simulations span 3.95 x 3.95 x 2.85 Mm 3 with a res- 
olution of 100 x 100 x 215. We started each simulation 



from the same state and evolved them until they achieved 
a new thermal relaxation based on their respective ra- 
diative models. Statistics were gathered after thermal 
relaxation. 

The following sections compare various space- and 
time-averaged quantities from the two simulations. Fig- 
ures show the mean values with standard error as a func- 
tion of height in the simulation domain. The zero point is 
defined as the location where the average temperature is 
equal to the effective temperature. We use the following 
definitions for statistical quantities. 

1 1 fto+T rL x pL y 

?(*) - \ \ \ qdydxdt, (9) 

1 lj X L/y J tg J Q J Q 

where the quantity has been averaged over a time T in 
a box of horizontal cross-sectional area L x x L y , and 
to is a time after the simulation has relaxed. The time 
average is over a sufficient interval to obtain statistical 
convergence. The RMS of the same quantity is: 

> q" = (7) - (q) 2 - (io) 

i 

[> 4.1. Radiative Efficiency 

The purpose of the radiative transfer method is to 
transfer energy within the simulation domain, and de- 
termine how energy is ultimately removed (the stellar 
flux) . The local contribution to radiative heating or cool- 
ing comes from applying the radiative transfer solver to 
the inhomogeneous medium. Different radiative transfer 
methods can remove the same total energy from a given 
stratification, but the energy may be removed in different 
ways. In this case, the effective temperature (measured 
as F cx cT^) remains unchanged, but the thermal struc- 
ture is different. This effect would appear as a different 
transition from convective to radiative energy transport 
and would affect the shape of the superadiabatic gradient 
(V — V a d). We refer to this as the radiative efficiency. 
The radiative efficiency determines the rate of transition 
from convective to radiative energy transport. A more 
efficient method will start to be super-adiabatic deeper 
in the star, and will have a higher peak in the superadi- 
abatic gradient. 

Figure [4] compares the radiative heating averaged over 
space (horizontal) and time. The integral of this func- 
tion over height is the total energy flux. Both methods 
have essentially the same energy flux (their effective tem- 
peratures differ by only 30 K, or 0.5%). The radiative 
transfer method has no effect on the adiabatic structure 
where the diffusion approximation is valid, so the 3D Ed- 
dington and ray integration methods have the same deep 
stratification. Convection carries the energy flux in the 
optically thick region, and is not directly affected by the 
differences in radiative transfer schemes. 

The superadiabatic gradient in Figure [5] shows a cor- 
responding change. A rising parcel of hot gas becomes 
superadiabatic at greater depth below the SAL in the 
ray method compared with 3D Eddington. At the peak 
of the SAL, however, the 3D Eddington method radiates 
more efficiently resulting in a higher peak. This means 
that below the SAL, the ray method is more efficient 
at transferring energy radiatively, but becomes slightly 
less efficient throughout and above the SAL peak. The 
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Fig. 4. — Radiative cooling as a function of simulation depth. 
The total radiative fluxes (the integral over height) are the same 
but the shapes of the curves are different. Lines show the space- 
and time-averaged mean, with bars showing the standard error. 
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Fig. 5. — Superadiabatic gradient as a function of simulation 
depth. The radiative transfer method affects the detailed shape of 
the superadiabatic gradients. 

differences are subtle but induce changes to some flow 
properties (Figure [8J and the mean stratification (Fig- 
ure [6]) in the superadiabatic part of the domain. The 
superadiabatic gradient of the 3D Eddington simulation 
is more negative throughout the atmosphere. 

4.2. Mean Stratification 

Figure [5] shows the density stratification as a function 
of temperature. Both simulations have the same strat- 
ification in the convective region and differ only in the 
SAL and optically thin atmosphere. The stratifications 
agree in the deep region because the two radiative trans- 
fer methods compute the same total radiative flux for a 
given stratification. The energy transfer in the deep re- 
gion is entirely convective and is indirectly affected by 
the radiative energy transport closer to the surface. The 
details of the radiative energy transport schemes are dif- 
ferent in each simulation, so the shallow region near the 
top has relaxed to a different state. The result is that the 
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Fig. 6. — The thermal structure (mean density vs temperature) 
for simulations computed using the 3D Eddington approximation 
and the ray integration method. The structures are the same in 
the adiabatic region but differ through the SAL and atmosphere. 

temperature and density structures of the simulations are 
distinct. 

The density and temperature profiles diverge as the 
temperature gradient becomes significantly superadia- 
batic. The densities differ by a factor of approximately 
1.2 within the SAL (at T/T cS w 1.3), and by a factor of 
approximately 1.5 in the atmosphere (at T/T c g ps 0.84). 
Note that this discrepancy is similar to the ra nge ex- 
hibite d in the comparison of solar simulations by iKupkal 
(2009) . The differences are likely caused by the different 
radiative efficiencies. While the two simulations have the 
same total energy flux (the integral of Q ra d over height), 
their cooling profiles in Figure 0] distinct. 

Simulations contain many numerical and physical as- 
pects that can affect the results. The thermal structures 
for solar simulations computed using different codes with 
a variety of physical inputs and parameters are indeed 
different. The stratifications in Figure [6] are within the 
range of what is seen when comparing res ults from diffe r- 
ent codes [for example see Figure 1 from iKup ka (2009)]. 
Except for the radiative transfer, all parameters that de- 
fine the two simulations are the same, so the change in 
stratification seen here is due entirely to the different ra- 
diative transfer schemes. Both radiative transfer meth- 
ods result in SAL structures that are realistic. Although 
we cannot determine which is more accurate, these re- 
sults provide a measure of the accuracy of 3D radiative 
transfer schemes in simulations. 

The thermal structure of the atmosphere is often pre- 
sented as the variation of temperature with optical depth. 
In stellar models, this so-called T-t relation is imposed 
as a boundary condition. Typically, the Eddington T-t 
relation is used for computing stellar models, i.e., 

T 4 = ^ ff (r + £), (11) 

which is purely radiative, and does not include the ef- 
fects of convective overshoot. This is one of the reasons 
that stellar models are inaccurate near the surface. Tak- 
ing T-t relations from simulations and using them in 
stellar models is a simple way to include overshoot. 

We find that the simulated T-t relation is sensitive 
to the radiative transfer scheme. While similar to Fig- 
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Fig. 7. — T-t relations for simulations computed with different 
radiative transfer methods. The atmospheric structures are very 
similar near r = 1, but show significant differences at low r. 

ure [6] which shows structural differences in the SAL, 
Figure [7] highlights the differences at small r. Sev- 
eral analytic and semi-empirical T-t relations are in- 
cluded in Figure [7] for reference. The simulations de- 
viate from the purely radiative Eddington T-r relation 
because they include the effects of convective overshoot 
and turbulence. These differences arc expected for realis- 
tic atmospheres , as de mons trated by the semi-em pirical 
iKrishna Swarnyl (|1966h and IVernazza et all (119811 ) mod- 
els. Simulations computed by iTrampedach eTin (fl997Tl 
also show changes in T-t relations as a result of convec- 
tion. 

One of the consequences of convective overshoot is a 
shift in the location of the mean T e ff surface with respect 
to mean optical depth. The T e ff surface from the 3D 
Eddington simulation is very close to r = 2/3, while the 
ray integration simulation is at slightly smaller (r s=s 0.6) 
optical depth. The slight differences in the turbulent 
properties of the simulations (which is a result of the 
different entropy profiles from using different radiative 
transfer models) produces slightly different positions of 
the T c ff surface. A shift in the T c g surface is not entirely 
surprising considering that semi-empirical T-t relations 
have mean T e ff surface at much smaller r. 

The simulations both produce plausible T-t relations. 
While distinct from each other, they are both within the 
range of solar semi-empirical T-t relations. Both simu- 
lations are closest to the analytic Eddington relation at 
optical depths near to or greater than unity, but deviate 
from it higher in the atmosphere. At small optical depth 
(r < 0.2) the r a y sim ulation is closest to the relation of 
IVernazza et al.l (|1981l ). while the 3D Eddington sim ula- 
tion appears to be similar to IKrishna Swamvl (jT966) . 

4.3. Turbulent Pressure and Gas Dynamics 

The convective velocities in the SAL result in added 
support in the form of turbulent pressure. Stellar models 
computed with mixing length theory do not include the 
effects of turbulent pressure. The added pressure from 
turbulence modifies the hydrostatic balance and is an 
important consideration for improving stellar models. 

Velocity statistics remain largely unchanged between 
the simulations. Figure [5] shows that the RMS of the 
vertical velocity as a function of height is nearly the 
same in both the ray and 3D Eddington simulations. The 
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Fig. 8. — The RMS of the vertical velocity remains essentially 
unchanged. Changing the radiative transfer scheme introduced 
only very small differences above the SAL. The surface is defined 
as height at which < T >=< T c g >, and the peak of the SAL is 
approximately 80 km below the surface. 

RMS vertical velocities are identical below the peak of 
the SAL, and show very minor differences above. This is 
because the velocities are generated in the adiabatic con- 
vectively unstable region, and both simulations have the 
same adiabat. As seen in Figure [SJ they can only start 
to deviate from each other in the SAL and convective 
overshooting regions where the different radiative trans- 
fer schemes have altered the stratification. 

Convection provides a significant source of support in 
the SAL. Turbulent pressure is defined as: 

fturb = pv% ms , (12) 

where uVms is the RMS of the vertical velocity. Turbu- 
lent pressure is small in the adiabatic region where gas 
pressure is high, reaches a maximum in the SAL, and 
drops quickly in the atmosphere where gas density is low. 
The two simulations show a relatively small change in 
-Pturb/^gas from 13% to 14% when the radiative transfer 
model is changed. Again, this is primarily a result of the 
differences in thermal structure (the density is different 
for a given temperature) since we see little change in the 
RMS velocity. 

While the radiative transfer scheme is the sole cause 
for the differences in Figure El turbulent pressure is also 
sensitive to other aspects of the simulation. As was the 
case with the mean stratification in Figure |H1 the change 
in turbulent pressure caused by the different radiative 
transfer schemes is within the range of what is seen in 
simul ations computed u sing different codes. For exam- 
ple, iBeeck et al.l (|2012t ) shows Pturb/-Pgas ranging from 
16% to 18% between solar simulations computed with 
the STAGGER, C05BOLD, and MURaM codes. 

5. CONCLUSIONS 

We used 3D radiation hydrodynamic simulations to 
compute the thermal structure near the surface of stars 
where energy transport transitions from fully convective 
to fully radiative. Since convection is multidimensional, 
achieving proper domain coverage is essential in provid- 
ing accurate radiative transfer. Radiative transport in 
3D is computationally demanding, and to keep compu- 
tational costs manageable, we make approximations to 
either the transfer equation or to domain coverage. 
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Fig. 9. — Turbulent Pressure as a function of height. Changing 
the radiative transfer scheme results in a small change in the peak 
turbulent pressure. Changes in turbulent pressure are primarily a 
result of the different mean stratifications, and not changes in the 
gas dynamics. The surface is defined as height at which < T >=< 
T c (f > , and the peak of the SAL is approximately 80 km below the 
surface. 

Our convection code includes two radiative trans- 
fer schemes that fundamentally differ in their meth- 
ods for domain coverage. The commonly used long- 
characteristic ray integration provides accurate solutions 
along individual rays, but lack of computational power 
necessitates under-sampling the domain and precludes 
proper 3D coverage. The highly variable and asymmetric 
nature of Q ra d makes it difficult to accurately sample the 
full 3D domain unless a great many rays are used, which 
is currently computationally cost prohibitive. Our expe- 
rience with undersampled ray integration suggests that 
it is important to consider the placement of the rays with 
respect to the highly asymmetric local flux, as the total 
flux may be sensitive to how well the local flux is sam- 
pled by the rays. We contrast the ray method with the 
3D Eddington approximation, which is derived from mo- 
ment equations closed at the level of first order, and has 
some computational advantages over long characteristic 
ray integration. 

We find good general agreement in mean quantities 
from simulations computed with the 3D Eddington ap- 
proximation and the long-characteristic ray integration 
method. Both simulations have the same adiabatic struc- 
ture in the convective interior, and the same energy flux. 
They do, however, differ in radiative efficiency as a func- 
tion of depth. The differences in energy transport pri- 
marily affects the thermal structure of the SAL and at- 
mosphere, while leaving the gas dynamics largely un- 
changed. 

The most significant differences we find are in the mean 
temperature stratification through the SAL and the T-t 
relation in the atmosphere. Although the T-t relations 
are distinct, they both compare well with semi-empirical 
solar T-t relations. Limb-darkening measurements may 
lead to improved constraints on the solar T-t relation. 
We note that results from ray integration and the 3D 
Eddington approximation may differ in the photosphere 
beyond the domain of the simulations studied here. Fur- 
ther study is necessary to determine whether there are 
differences in the photosphere of the 3D RHD simula- 
tions, and to what extent any differences have on ob- 
servable quantities, such as features derived from spec- 
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tral line synthesis using structures computed from the 
3D simulations. 

Simulations computed with the 3D Eddington approx- 
imation and under-sampled long-characteristic ray inte- 
gration both provide similar thermal structures that self- 
consistently couple the region of efficient adiabatic con- 
vection to the atmosphere. Indeed, the differences seen 
by changing the radiative transfer scheme are within the 
range of resu lts of solar simulations compute d using dif- 
ferent codes (|Kupkall2009t iBeeck et al.ll2012f ). We there- 
fore conclude that both radiative transfer techniques give 
rise to robust results for the angular resolution afforded 



by current computational power. 
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